First, use this code to download the data from the USGS site.

library(dataRetrieval)
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(tidymodels)
── Attaching packages ────────────────────────────────────── tidymodels 1.3.0 ──
✔ broom        1.0.7     ✔ rsample      1.2.1
✔ dials        1.4.0     ✔ tune         1.3.0
✔ infer        1.0.7     ✔ workflows    1.2.0
✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
✔ parsnip      1.3.1     ✔ yardstick    1.3.2
✔ recipes      1.2.0     
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter()   masks stats::filter()
✖ recipes::fixed()  masks stringr::fixed()
✖ dplyr::lag()      masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step()   masks stats::step()
library(tsibble)
Registered S3 method overwritten by 'tsibble':
  method               from 
  as_tibble.grouped_df dplyr

Attaching package: 'tsibble'

The following object is masked from 'package:lubridate':

    interval

The following objects are masked from 'package:base':

    intersect, setdiff, union
library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout
library(feasts) #autoplot visualizations
Loading required package: fabletools

Attaching package: 'fabletools'

The following object is masked from 'package:yardstick':

    accuracy

The following object is masked from 'package:parsnip':

    null_model

The following objects are masked from 'package:infer':

    generate, hypothesize
library(fable)

#If you prefer to keep the date as a date object, you can use lubridate::year() and lubridate::month() to extract the year and month, and then combine them:



poudre_flow <- readNWISdv(siteNumber = "06752260",    # Download data from USGS for site 06752260
                          parameterCd = "00060",      # Parameter code 00060 = discharge in cfs
                          startDate = "2013-01-01",   # Set the start date
                          endDate = "2023-12-31") |>  # Set the end date
  renameNWISColumns() |>                              # Rename columns to standard names (e.g., "Flow", "Date")
  mutate(Date = as.Date(Date, format = "%Y-%m")) |>  # Create a Year-Month column as "YYYY-MM"
  group_by(Date) |>                         # Group by the new Year-Month column
  summarise(Flow = mean(Flow, na.rm = TRUE))           # Calculate the average daily flow for each month
GET:https://waterservices.usgs.gov/nwis/dv/?site=06752260&format=waterml%2C1.1&ParameterCd=00060&StatCd=00003&startDT=2013-01-01&endDT=2023-12-31

Convert to tsibble: Use as_tsibble() to convert the data.frame into a tsibble object. This will allow you to use the feast functions for time series analysis.

# Assuming you already have 'poudre_flow' with the "Date_YearMonth" column as shown earlier
poudre_flow_tsibble <- poudre_flow %>%
  as_tsibble()  # Convert to tsibble using Date as the time index
Using `Date` as index variable.

Plotting the time series: Use ggplot to plot the time series data. Animate this plot with plotly

p_plot<- poudre_flow_tsibble %>%
  autoplot() +
  geom_line(color= "purple") +
  labs(
    title= "Interactive time series plot of the Poudre Rive flow",
    x= "Date",
    y= "Flow")
Plot variable not specified, automatically selected `.vars = Flow`
ggplotly(p_plot)

Decompose: Use the model(STL(…)) pattern to decompose the time series data into its components: trend, seasonality, and residuals. Chose a window that you feel is most appropriate to this data. Describe what you see in the plot. How do the components change over time? What do you think the trend and seasonal components represent?

# STL decomposition model doesn't assume constant seasonal effect like decompose()

# Convert to a regular time series object for STL (monthly data)
flow_ts <- ts(poudre_flow_tsibble$Flow, start = c(2013, 1), frequency = 12)

# Apply STL decomposition
stl_decomposed <- stl(flow_ts, s.window = "periodic") %>%
  plot()

This plot shows the seasonal trends from 2013 and predicts them up until about 2350. They follow typical oscillations that occur on a yearly basis, except they slowly start having lower maximum flows. The trend could represent the effects of climate change on the river flows over time, and the seasonal component is driven by the typical natural biological processes as well as climate change (some seasons may start to get slightly shorter).